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Abstract: Discrete fine-scale models, in the form of either particle or lattice models, have been formulated suc¬ 
cessfully to simulate the behavior of quasi-brittle materials whose mechanical behavior is inherently connected to 
fracture processes occurring in the internal heterogeneous structure. These models tend to be intensive from the 
computational point of view as they adopt an “a priori” discretization anchored to the major material hetero¬ 
geneities (e.g. grains in particulate materials and aggregate pieces in cementitious composites) and this hampers 
their use in the numerical simulations of large systems. In this work, this problem is addressed by formulating a 
general multiple scale computational framework based on classical asymptotic analysis and that (1) is applicable 
to any discrete model with rotational degrees of freedom; and (2) gives rise to an equivalent Cosserat contin¬ 
uum. The developed theory is applied to the upscaling of the Lattice Discrete Particle Model (LDPM), a recently 
formulated discrete model for concrete and other quasi-brittle materials, and the properties of the homogenized 
model are analyzed thoroughly in both the elastic and inelastic regime. The analysis shows that the homogenized 
micropolar elastic properties are size-dependent, and they are functions of the RVE size and the size of the ma¬ 
terial heterogeneity. Furthermore, the analysis of the homogenized inelastic behavior highlights issues associated 
with the homogenization of hue-scale models featuring strain-softening and the related damage localization. Fi¬ 
nally, nonlinear simulations of the RVE behavior subject to curvature components causing bending and torsional 
effects demonstrates, contrarily to typical Cosserat formulations, a signihcant coupling between the homogenized 
stress-strain and couple-curvature constitutive equations. 
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1 Introduction 


Discrete fine-scale models, in the form of either particle or lattice models, have been formulated successfully 
in the literature to simulate the behavior of a variety of different materials. Their use has become more and 
more popular in the last few decades due to a number of appealing properties that make them advantageous 
compared to continuum based formulations. 

The geometry of discrete models is built with reference to the actual internal structure of the material 
of interest and it consists of “particles” connected through either “contact points” or “connecting struts” 
(also called “lattice elements”). This “a priori” discretization allows simulating material heterogeneity 
efficiently in the case of materials - such as concrete, rock, sea-ice, and toughened ceramics - characterized 
by hard and stiff inclusions embedded in a more compliant, weak, and brittle, matrix. In addition, the 
intrinsic particle/lattice spacing automatically provides the formulation with an internal characteristic 
length which can be made randomly variable if the discrete model is constructed according to the actual 
random distribution of material heterogeneity. 

The degrees of freedom (displacements and rotations) are defined only at a finite number of points - 
referred also as “nodes” thereinafter - which, depending on the formulation, may or may not correspond 
to the partice center of mass or particle centroid. Strain and stress measures are defined at a finite 
number of points coinciding with the contact points or with some specified points along the connecting 
struts. The constitutive behavior is formulated through vectorial, as opposed to tensorial, stress versus 
strain relationships and stress tractions are supposed to be distributed over either a “contact area” or 
the cross sectional area of the connecting struts (in this paper, this area will be generically referred 
to as “facet”). Finally, the classical concepts of equilibrium and compatibility are formulated through 
algebraic equations, instead of partial differential equations typical of continuum mechanics. One of the 
main advantages of discrete models is that the discreteness of the formulation permits handling naturally 
displacement discontinuities arising during damage localization and fracture processes. 

Rigid particle models, under the name of Discrete Element Method (DEM), were first formulated to 
simulate both natural materials, such as geomaterials p El El m, as well as man-made materials like 
concrete mum- A somewhat similar model is the rigid-body-spring model (RBSM), which subdivides 
the material domain into rigid polyhedral elements interconnected by zero-size springs [lElISlIin!. 
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Lattice models, pioneered by Hrennikoff m to solve elastic problems in the pre-computers era, were 
later developed by many authors to model fracture in quasi-brittle materials in both 2D [12], and 3D 

[Mlliaiisillel. 

More recently, various discrete models, in the form of either lattice or particle models, have been quite 
successful in simulating concrete materials Ha [HI HH na ED]. For an extensive review of the currently 
available models for concrete the reader is directed to a recent special issue [2T] collecting several papers 
covering a wide variety of concrete mechanics phenomena spanning several length scales, from the scale of 
cement particles to that of reinforced concrete structural members. 

In most applications of interest in practice, fine-scale models lead to fairly large computational systems 
characterized by a huge computational cost making their practical use rather limited. For example, the 
full-scale computational analysis of an average concrete bridge would require millions of degrees of freedom 
or the simulation of a rock formation would require billions of degrees of freedom. The solution of such large 
problems, although possible in principle with large super computer clusters, is unimaginable in everyday 
engineering practice. For this reason, many studies have been devoted to finding optimal and rigorous 
approaches for multiscale computation. 

Among different multiscale techniques available in the literature [22] , the ones based on homogenization 
theory have been widely used over the past decades. The homogenization theory relies on two main 
assumptions. The first is the existence of a certain volume of material, the so called Representative Volume 
Element (RVE) or Unit Cell (UC), carrying a complete description of the internal material structure 
[2S1I21]. The second is that the size of such a volume is much smaller than the size of the overall solid 
volume under consideration. The latter is also known as the “scale separation” assumption. 

Hill [2S] , Eshelby [2S] , Hashin and Strikman [27] pioneered analytical homogenization techniques which 
were developed later by other authors [2HI ES] • Analytical homogenization is able to reasonably approxi¬ 
mate material properties when the exact solution of the boundary value problem associated with the RVE 
problem can be obtained. However, in this approach, elastic behavior, small strains, and relatively sim¬ 
ple internal structure are the limiting assumptions typically adopted. When complicated heterogeneous 
structures are considered, or constitutive behavior of constituents are nonlinear, other homogenization 
techniques [301 El] needs to be considered. 


3 


To overcome these difficulties, computational homogenization is often used in the literature |321 EHl 
EH EH- III fhis approach, a single RVE is assigned to each calculation point (e.g. gauss point in a Finite 
Element mesh) in the macro domain and at each step of the nonlinear analysis, macro-strain increments are 
imposed as essential boundary conditions to the RVE. The solution of the RVE boundary value problem is 
then averaged for the calculation of the associated macroscopic stress tensor. Since no assumption is made 
for the macroscopic constitutive law, this method can be used for materials featuring extremely nonlinear 
behavior. 

A somewhat similar but more mathematically rigorous homogenization technique is the so-called 
Asymptotic Expansion Homogenization (AEH) that uses the asymptotic expansion of the displacement 
field based on a length parameter representing the ratio between the length scale of material heterogeneity 
and the macroscopic length scale. Starting from this expansion hierarchical boundary value problems are 
obtained at different scales. This approach can easily handle problems with multiple (more than 2) scales 
in both space and time [35] ; it does not make assumptions on the character of the macroscopic constitutive 
equations; and its implementation in computer codes is relatively simple. 

Within the extensive literature on AEH, remarkable is the work of the following authors. Hassani 
[3611371 investigated formulation of homogenization theory and topology optimization and its numerical 
application to materials with periodic microstructure. Chung [SH| presented detailed derivation of multiple 
scale formulation for elastic solids. Fish employed this approach to study elastic as well as elasto-plastic 
composites [39] . Ghosh [30| adopted MH along with Voronoi Cell Finite Element Method (VCFEM) to 
study the behavior of composites with random meso-structure m More recently, Fish [HS] introduced 
the Generalized Mathematical Homogenization (GMH) to derive continuum constitutive equations starting 
from Molecular Dynamics (MD). 

All the aforementioned work is relevant to Cauchy continuum formulations. However, homogenization 
schemes were also used for the multiscale analysis of Cosserat continuum models, in which an independent 
rotation field appears in addition to the displacement field. Feyel [HS] built a homogenization scheme to 
couple a Cauchy continuum formulation at the micro-scale giving rise to a Cosserat continuum formulation 
at the macro-scale. Asymptotic homogenization technique was employed by Forest [12] for upscaling 
elastic Cosserat solids. In this work, the author studied various types of asymptotic expansions for the 
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displacement and rotation fields and investigated their effect on the resulting macroscopic continuum 
behavior. Results of this investigation, showed that the nature of the homogenized continuum depends on 
the ratio of the Cosserat characteristic length of constituents, size of heterogeneity and typical size of the 
structure. 

Chan et al. jl3] derived the governing constitutive equations for strain gradient elasticity for both 
homogeneous and functionally graded materials using the strain energy density function and the related 
definitions of the stress fields. They showed that additional terms appear in the equations that are related 
to the strain gradient nonlocality and the interaction between material nonhomogeneity. Bardenhagen et 
al. [Hj obtained a nonlinear higher order gradient continuum representation of discrete periodic micro¬ 
structures by means of an energy approach. The developed model was then employed to investigate the 
existence and stability of localization bands and their relationship to the model loss of ellipticity. Finally, 
homogenization of discrete atomic models into equivalent continuum can be found in publications where 
the authors exploited asymptotic analysis techniques na and the mathematical F-convergence method 

m- 

The present study derives a general multiscale homogenization scheme suitable for upscaling materials 
whose fine-scale behavior can be successfully approximated through the use of discrete models featuring 
both translational and rotational degrees of freedom. 

2 The Fine-Scale Problem 

With reference to Figure [^, let us consider the interaction of two adjacent particles, I and J, sharing 
a generic facet. If one limits the analysis to the case of small strains and displacements - which is a 
reasonable assumption in absence of large plastic deformation prior to fracture as observed in brittle and 
quasi-brittle materials - meaningful measures of deformation m can be defined as 

= 1 (U-^ + X c-^ - X c^) • (1) 

and 

x" = ; (e" - ®') ■ e" (2) 
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where = facet strains; = facet curvatures; r = — x-^ is the vector connecting the 

particle nodes P/ and Pj; {a = N, M, L) are unit vectors defining a facet Cartesian system of reference 
such that ejy = is orthogonal to the facet and ej/ ■ x^'^ > 0; U^, = displacement vectors of node P/ 

and Pj] 0^, 0“^ = rotation vectors of node Pi and Pj; and = vectors connecting nodes Pi and Pj 

to the facet centroid, see Fig. It must be observed here that displacements and rotations are assumed 
to be independent variables. 

For given strain and curvature vectors, a vectorial constitutive equation provide stress, and couple, 
tractions on each facet. Formally one can write = ta{eN, = ma{xN, where, 

in general, summation rule applies over a. As an example, the elastic behavior can be formulated through 
the following equations 


A^, Af, Z/^ 


(3) 


in which each traction component is proportional to the associated strain or curvature (summation rule 
does not apply); and Ba, Wa are fine-scale elastic constants which are related by a characteristic length 
£a- An example of nonlinear facet constitutive equations is reported in Appendix]^ Section p\..3[ with 
reference to the so-called Lattice Discrete Particle Model (LDPM) that will be considered in the numerical 
examples. 

Finally, the computational discrete fine-scale framework is completed by imposing the equilibrium of 
each single particle subject to the effect of all surrounding particles. Translational and rotational dynamic 
equilibrium equations read 


a/'u' + m;,0 -Vb" = y]>it 


ij 


(4) 


and 


M^0^ = + m 


IJ\ 


(5) 


where x is the moment of the traction with respect to the particle node Pi; Pi is 

the set of facets surrounding node P/ and obtained by collecting all the facets associated with each node 
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Figure 1: Geometrical explanation of the two-scale problem: (a) Geometry of two neighboring particles, (b) 
Macro material domain, (c) Meso-scale domain with material heterogeneity. 


pair (/, J); A = facet area; superimposed dots represent time derivatives; is the particle volume; b° 
is the body force vector; = mass of node Pj] and = moment of inertia tensors. It is worth 

observing that = 0 and = M/l if the particle node is the particle center of mass; the axes of 
the system of reference are parallel to the particle principal axes of inertia; and the principal moments 
of inertia are the same in all directions. These conditions, although applicable only to a limited number 
of cases (e.g. spherical particles), do not reduce the conceptual generality of the derivation that will be 
presented in this paper and will be assumed thereinafter for simplicity. 

3 Asymptotic Expansion Homogenization 

In this section, the two-scale homogenization of the general fine-scale problem introduced in the previous 
section is pursued by means of the approach proposed in Ref. [3H] • In the original formulation only central 
forces were assumed to act on the particles and, consequently, the rotational equilibrium equation was not 
considered. 

3.1 Two Scale Approximation and Asymptotic Expansions 

In order to perform a two-scale asymptotic expansion homogenization, a periodic discrete system, composed 
by a number of adjacent RVEs, is considered in this section. In Figure[l|3, the generic macroscopic material 
domain and the corresponding global coordinate system X are shown. At any point in the macroscopic 
domain, two separate length scales and the corresponding local coordinate systems, x and y, are introduced 
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to represent ( 1 ) the macroscopic domain, in which the problem is defined as homogeneous continuum with 
no detail of material heterogeneity, and ( 2 ) the meso-scale domain, in which heterogeneity is modeled by the 
discrete meso-scale model. Vector X, as shown in the figure, is the vector connecting the origin of the global 
macroscopic coordinate system to the mass center of a generic RVE. In Figure [^, a zoomed view of the 
macroscopic material point is shown in the local meso-scale coordinate system y, in which a representative 
volume of heterogeneous material is depicted. One should consider that in Figure [T^, particles / and J 
are shown in the local macroscopic coordinate system x. Therefore, they should be plotted in smaller size 
compared to Figure [T]:, but this was not done for the sake of clarity. If the separation of scales exists, one 
can write the following relationship linking macro and meso local coordinate systems 

X = rjy] 0 < 77 << 1 ( 6 ) 

where 77 is a very small positive scalar. In addition, the displacement of a generic node P/, = 

u(x^,y^), can be approximated by means of the following asymptotic expansion 

u(x,y) u°(x,y)+ 77u^(x,y) (7) 

where only terms up to order 0 ( 77 ) are considered. Functions u°(x, y), and u^(x, y) are continuous 
with respect to x and discrete (i.e. defined only at finite number of points) with respect to y. 

In order to define the asymptotic expansion for rotations, it is convenient first to postulate the existence 
of a continuous displacement-like field d’^(x) such that 20-^ = V x d’^|x=x^- If d’^(x) is replaced by a two- 
scale approximation similar to the one in Equation]^ one can write, 0^ = 0(x^,y^), and 

6>(x,y) ^ 77-^a;°(x,y) + (^°(x, y) + a;^(x, y) + 77(^i(x,y) (8) 

where 2u;° = Vj^ x d°; 2y)° = x d°; 2u?^ = Vj^ x d^; 2y)^ = x d^; and subscripts x and y 
identify the nabla operator in the coarse- and fine-scale, respectively. Thus, u;°, should be interpreted 
as rotations in the fine-scale whereas as the corresponding coarse-scale rotations. It is worth 

observing here that, contrarily to the expansion of displacements, the asymptotic expansion for rotations 
features a term of order 0 ( 77 “^) and two distinct terms of order (9(1). 


In the macroscopic coordinate x, the difference in position between nodes Pj and Pj can be considered 
as infinitesimal. Hence, in order to obtain the asymptotic expansion of strains and curvatures, it is 
convenient first to obtain the Taylor series expansion of displacement and rotation at nodes Pj around 
point Pj of coordinate in the local coordinate system x. By assuming that the displacement and rotation 
fields in Equations and [8| are continuous and differentiable with respect to x, one can write 


U- = «i(x-', y") = ui + uP xY + xYxY + • • • 


(9) 


©i =^i(x ,y ) = 0 ^ + 9 ijXj +- 9 ijkXj x^ + 


( 10 ) 


where uf = Ui(x^,y^); uP = dui/dxjipP = d^Ui/dxjdxk{:>^^ ,y^)] 6*/ = 6'j(x^,y^); 9P = 
d9i/dxj(yJ9f-j^ = d^9i/dxjdxk{'^^,7'^)', xY is a vector connecting node Pj to node Pj in the x space. 
By substituting Equations and into Equation [T} and using the Taylor expansion of displacement and 
rotation of node Pj around node Pj (Equations and 10) one obtains the multiple scale definition of facet 
strains (see Appendix for details) 


^a = ri + e° 

where 


(11) 


e-' = f-' 


Ui Ui + SijkiOj 


JJ 


( 12 ) 


e“ = f-‘ 


u. 


Ij I ..1/ 


, 0 J . IJ . , 0 J \ -J 


i + or - < + w + ^r + - ^^3k 


UI I , , 1 / \ -I 


JJ 


(13) 


= r ^ 


n,JnJJ I r, 0 J I I r JJ _L r JJ n J J _L / ^J o J J _L ^ , 0 ^^ nJJ^JJ \ J r J^ J 

^ijyj + O^iJkyj yk + ^ijk l Jj + Jj,mym + ^j,mym + Tl^j.mnym yn ] ^k ~ ^ijkjj 


JJ 


(14) 


In the previous equations, Sijk is the Levi-Civita permutation symbol and length type variables have 
been changed into their 0(1) counterparts by using Equation]^ r = rjf, = rjc'l. 

Similarly, multiple scale definition of facet curvature can be calculated as (see Appendix for details) 
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VXa = V + C + #« 


(15) 


where 


=r ^ 


, , ,0-f 


JJ 


€ = r-^ 


9 ?°-^ + uj]-’ + - cjI^ 


JJ 


J = r ^ 


.JJ _L / ^JnJJ _L .JJn.IJ _L . .^J n.IJn.IJ . JI 

+<^i,jyj +^ijyj +2 i’^kyj yu -^i 


JJ 


(16) 


(17) 


(18) 


It is worth noting that in this section as well as in the rest of the paper superscript IJ has been dropped 
when the permutation of / and J is not associated with a sign change. 


3.2 Multiple-Scale Equilibrium Equations 

In order to obtain the correct scale separation of the governing equations, a rescaling of the discrete 
equilibrium equations needs to be performed. For the sake of simplicity, and since only quasi-static 
problems are concerned in the current research, it is assumed = 0 and = M/l on the left hand 
side of Equation |4| Rescaling is pursued by assuming that the material density, mass per unit volume, 
is of order zero: p ~ 0{1), which along with the displacement asymptotic expansion implies that the 
left-hand-side of Equation]^ is ~ 0{ri^). By dividing both sides of Equation]^ by and considering that 
all length variables should be considered ~ one obtains 


Min’ -V‘h'> = r,-'J2At^e" (19) 

where = M^/rf, /rf , A = A/rf are all quantities ~ 0{1)- For reason of dimensionality, 

body forces 6 ° can be always assumed to be proportional to gravity pg and, consequently, they can be 
considered 0(1) quantities as well. 

One can rescale the rotational equation in a similar fashion by recognizing that, according to the 
previous discussion, the rotational moment of inertia is ~ 0(ri^). Dividing both sides of Equation]^ by 77 ^ 
one obtains 
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( 20 ) 


riMfe^ = ri ^^A{ri ^ujae^^ + r] 

J^i 

where M/ = Mq /rf is ~ 0{1). 

In the elastic regime one can write: ta = V~^ta^ + + Vta'^ where ta^ = Eaea\ and Ea is assumed to 

be ~ 0(1). In addition, + ?q + ^?q which gi'^ = Wa'ipa'^] Wa = Ea^a'i ^ 

Finally, Pa = V~^'^a = V~^Pa^ +Pa + VPa where Pa^elf = x . Since Wa and are moments, it is 

reasonable that the asymptotic expansion of those variables divided by r] is similar to the one for tractions 
ta, considering that length type variables are considered to be r\j 0(r,). 

Introducing these traction expressions along with the asymptotic definition of displacement and rotation 
fields Equations and (which also imply rjd^^ = + 0{p)), into the rescaled equilibrium equations 

leads to 


-2 


Ti 


+ p 


Ti 




Tl 


JJ 


+ 0(p) - 0 




( 21 ) 


and 


:Fi Ti 

+ 5] 24 (p^^ + g^e^") + 0(p) = 0 

Ti 

in which terms of different orders are gathered together. The multiple scale equations reported above 
can also be used for nonlinear constitutive equations provided that facet tractions and facet moments can 
be expressed through the multiple scale decomposition exploited above. It will be shown later in the paper 
that this can be indeed achieved under some reasonable assumptions. 



3.3 The RVE Problem 


Let’s first consider the equilibrium equations at the 0{p scale. From Equations 21 and 22, it is evident 
that the 0{p~‘^) equilibrium equations represent the equilibrium of all particles in the RVE subjected to the 
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stress tractions ^ and the moment tractions ^ and without any applied external load. Consequently, 
solution of the 0{ri~^) problem implies t~^ = 0 and q~^ = 0 , which in turn, leads to = 0 and 'tp~^ = 0 . 


By taking into consideration the dehnitions of and (Equations 12 and 16) such result indicates 
that the 0{ri~‘^) problem represents a rigid body rototranslation of the RVE. This can be expressed as 


u°(X, y) = u°(X) + £q,|/,a;°(X) (23) 

in which the helds and are only dependent on macroscopic coordinate system X, i.e. these 
quantities varies smoothly in the macro-scale material domain; they do not change within the RVE domain; 
and they can be calculated when kinematic boundary conditions are specified for the 0{ri~^) problem. 
These boundary conditions must describe the physical fact that the RVE is attached to a point in the 
macroscopic continuum. Hence, v° must correspond to the macroscopic displacement field, and must 
be equal to the macroscopic rotation field: Since is constant over the RVE, then is also 

constant in the RVE. 

On the basis of Equationand the discussion above, one can rewrite the 0(1) strains and curvatures 
as (See Appendix [C| for details) 


= r' {u'/ - ui‘ + - SijtuYy) e'y + fg (7,,+ 

</.“ = f-' (u^y - yy e'j+ p-k,, 


(24) 

(25) 


where 7 *^ = K,ij = are the macroscopic Cosserat strain and curvature tensors, respec¬ 

tively. The vector is the position vector of the centroid of the common facet between particle / and J 
and Pj" = is a projection operator. Comparing the first term of Equation 24 with Equation it 

can be concluded that this term is the lower scale definition of the three components of the facet strains 
(one normal and two tangential) written in terms of fine-scale displacements and rotations and a;^. 

-fq (7q Is the projection of macroscopic Cosserat strain 


The second term of Equation 


24 


and curvature tensors on each facet. Similarly, Equation 25 shows that the 0(1) curvature includes a 
fine-scale term (see Equation]^, which depends on fine-scale rotation term and a coarse-scale term 
corresponding to the projection of macroscopic curvature tensor on each facet. Therefore, Equations 1^ 
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and 25 express the 0(1) facet strains and curvatures as the sum of their fine-scale counterparts and the 
projection of macroscopic strain and curvature tensors onto the facet level. It is worth nothing that the 
projection operator corresponds exactly to the one used in the microplane model [IZlIlHj if 
i.e. the discrete model is formulated in such a way the facets are orthogonal to the associated lattice struts. 
In addition, it must be noted that the term SjmnK-imyn transforms the macroscopic curvature tensor, which 
is constant over the RVE, to different strain values at different positions inside the RVE, which is 
then projected on the facets through the operator P^. Expanding this term for different components of 
curvature tensor, it can be shown that it perfectly corresponds to the strain field generated by curvatures 
in classical beam theories. 

Strains and curvatures of order 0{ri) can also be rewritten by taking into account Equation One 
gets 


e' = f-^ 


uljvy + SijwYci + SijkU^YmyY4 - eijkY/ci 


1 


+ -TyvljkyYvk + -Ty^ijk^lmnyt;! ynyl + ^ijki^lmyt^ 4 


1 




,0 


JJ 


4a = r 


^-1 


^AJ I , YJ^JJ ,AI I , ,0 ^Jj I , ,0J nJJnJJ 


^OLl 


(26) 

(27) 


Detailed mathematical derivation of Equations through is provided in Appendix 
In the previous derivation, where linear elastic behavior was assumed, the equilibrium equations at the 
0{r]~^) scale were shown to represent the rigid body motion conditions for the RVE and, consequently, 
they led to zero strains, e“^, curvatures, tractions, and moments, and q^, at the 0{riA 
scale. These conditions can be reasonably assumed a priori in the case of nonlinear material behavior. 
In this case case one may write ta = Pa = Pa{4 + ^a = Qa{v~^4 which 

a, P = N, M, L. Since 77 is a small quantity, one can also write the Taylor expansion of ta and Pa around 
the 0(1) component of strain and the Taylor expansion of around the 0{r]~4 component of curvature: 
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ta = tai4 + ^ 4 ) = ^ 4 


Pa = Pai^^H + ^4) = P«(^S) + ^ 


^Pa(e°) 1 

aeO 


( 28 ) 


.-i./.o , . -i.,, 0 N , 


= qaiV + "0/?) = 5a(?? 0^) + P- 


5^0 


-0. 


which can be rewritten as ta = t° + ptl^; Pa = Pa + VPa] Qa = Qa~^ PQa^ with the following conditions 


= ^a(e?); Pa = Pa(e^); = ?«(?? V°); 

1 _ ,1 

'S’ ~ 


1 ^1 ^P°n. .1 


(29) 


= 

a Q 0 7’ 

S """T 


“ deO aeO 


This demonstrates that Equations and ^ are valid also in the case of nonlinear material behavior 
under the assumption that traction and moments at the 0{ri~^) scale are zero as required, in the linear 
case, by the rigid body motion of the RVE. 


The RVE problem is governed by the 0{ri ^) terms in Equations 21 and 22, Considering those terms 
and scaling back all the variables, one can write the 0{ri~^) equations as 


= Y:Mc'-ty/ + myj) = o 

Ti Ti 


(30) 


Equations are force and moment equilibrium equations of each single particle inside the RVE sub¬ 
jected to (9(1) facet traction and moment vectors, which, in turn, are functions of e° and 
consisting of a coarse-scale and a fine-scale term (see Equations [24| and 25). In other words, Equations 30 
state that the macroscopic strain, 7 jj = uV — and curvature, = a;V, tensors should be applied 

on all RVE facets as negative eigenstrains, and the fine-scale solution, in terms of displacements u\ and 
rotations u} of each particle, must be calculated satisfying its force and moment equilibrium equations, 
while periodic boundary conditions are enforced on the RVE. The solution of the equilibrium equations 
also provides facet traction and moment vectors that are later used to compute the macroscopic 
stress and couple tensors. 
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3.4 The Macroscopic Problem 


Finally, let us consider the 0(1) equilibrium equations in Equations 21 and 22 The 0(1) translational 
equilibrium equation for each particle in the RVE reads 








del 


e^ + 




(31) 


where all the variables have been scaled back in the original system of reference, and 
By using Equation and by averaging the contribution of all particles in the RVE, one can write (see 
Appendix for details) 


”0 ^ A 1 7 


(32) 


1 


where Vq is the volume of the RVE; = X/ -Vl/K) is the mass density of the macroscopic continuum; 
bi = 6°(1 — (/)); and 0=1 — '^jV^/Vq is the porosity of the macroscopic continuum. Equation 


32 


was 


derived under the assumption that X/ = 0^ which corresponds to the assumption that the local 

system of reference is the mass centroid of the particle system within the RVE. 

Before proceeding with the derivations, let’s take a closer look at the definition of and the term 


/del)el on the RHS of Equation 32, Each facet in the material domain is shared between two 
particles, say I and J. Therefore, by summing up the contributions of two adjacent particles, one obtains 


dtl-^ 1 dtf 1 
*-ei + ^ei=- * 


del 


del 


r del 


^nW/ + £njk^y4 + ^njkytyyK 


+ \vl,jkVyVk’ + ^^njkl^lmoVy y^Vk + ^njk^lmyy4 ) ^an 


1 dty 

H- — 

r del 


U^yf + £njk<fY4 + er,jkUj]yyyy - Snjky/4 


2^^njkyj Vk 2^^'bijk^j^rnoym yo yk ^njk^j^rnym ^k j^an 


(33) 


Considering the definition of the vector = Y’’ — ■, one can write = —y^ and cl — y = y^. 

In addition, and hold for each facet. Finally, the sign of e° does not change by 
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interchanging I and J in its definition. This leads to = —dtf^/de'^. Taking all above facts into 

account, Equation can be written as 


dtfJ 1 dtf 1 lat" 
* -ei + -^ei = - * 


de^ 


deo 


r de^ L 


cfm \ n.m 


Un,m + ^njk(^]%ci - 


+ - ^njkOjl^ + enjki^lmoyi'’Vk ) e^n 


(34) 


Comparing the expression inside the bracket on the RHS of Equation 34 to the definition of in 
Equation it can be concluded that 


+ ^,1 = 

~l~ rx n r\ n r\ Um O c/m 


de^ 


de^ 


del 


dXr 


(35) 


Therefore, one can average the term /del)el on each facet and replace it with l/2{dtl'^/dxm)y-, 
in the equilibrium Equations which can be rewritten as 


ij 

m 


PuVl 


=—yyAr 


I Ti 


dtj'^ jj 
-K — ^7 + h 

OXn 


(36) 


Einally, by considering that (1) 71^)/dxj = dtj'^/dxjny + tl'^dny/dxj and (2) dny/dxj = 0 for 

the periodicity of the problem; and by recalling that , one obtains 


PuVi = cr“ , + bi 


(37) 


and 


2vi, E E 


(38) 


/ 


Equation is the classical partial differential equation governing the equilibrium of continua whereas 
Equation provides the macroscopic stress tensor by averaging the solution of the RVE problem. It is 
worth mentioning that Equation coincides with the virial stress formula for atomistic systems derived 
in Ref. j35], but it is also equivalent to the averaging formula used in the classical microplane model 
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formulation and derived through an energetic equivalence. 

The 0{1) moment equilibrium equation is considered next. Since the purpose in this section is to 
average the equation of motion of all particles inside the RVE and derive the macroscopic equilibrium 
equation governing the entire RVE, to have a consistent formulation for all particles and RVEs, one must 
consider the moment of all forces with respect to a hxed point in space. 

Eor the generic particle I, by taking the moment of all forces with respect to the origin of a global 
macroscopic coordinate system as shown in Figure[^ and by considering the results of the problem, 

one can write (see Appendix for details) 


(ii» 


^kmnV 


■Pi ^ 


de^ 


dml'^ 

d'lpa 




(39) 


where Xj is the position vector of particle I in global coordinate system; = SijkX^tlf is the 
moment of facet traction with respect to the point O; is the position vector of the contact point C 
between the particles I and J in the global coordinate system, and ml'^ = Also, xj and are 

the position vectors of the particle I and the contact point C with respect to the mass center of the RVE, 
respectively. 

By summing up the moment equilibrium equations of all particles inside the RVE and dividing by the 
volume of the RVE, and considering that = Xj + xj, one obtains (see Appendix [d| for details) 




(40) 


where + M^eijkekmnXjxi] /ho is the inertia tensor of the RVE. In deriving Equation 

40, the particle density /V^ was assumed to be constant for all particles; and the local system of 
reference at the center of the RVE was chosen such that M/x/xj = 0 for any i ^ j, i.e. as mentioned 
earlier in this paper, the axes of the system of reference are principal axes of inertia for the system of 
particles within the RVE. 

Before moving forward with the derivation, let’s first consider the second term on the RHS of Equation 
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40 


For a facet in the material domain which is shared between particles / and J, by summing the 


contribution of two particles / and J on the term and by considering the definition of 'tpl, 

(see Equation , one gets 


dml'^ 




i’, 


dm'p , 1 1 dml'^ 

-^a = 




0 
a 


r dip\ 


^ijyj ^ijyj o^hjkyj yk 


AJn.IJ 


a/ 


0 ../J 


,0^ nJJn.IJ 


JJ 


1 dm{^ 
f d'lp^ 




+ o^^idky^Vk 


,0 „,J/ 


1 


,07 „,JI„,JI 


(41) 




Since the moment stress vector applied on a single facet belonging to two particles / and J are the 
same in magnitude but opposite in direction, one can write mj'^ = —mf^; and consequently, dnij'^/'ip^ = 
—dmf^/'ip^. In addition, the sign of does not change by interchanging / and J in its definition, and 
that 


,iJ = 

im 


'Vm ^ ^cd — ~^od^ Equation 41 can be written as 


dmp'^ , 1 dmf^ ,, 1 dmp'^ 


dxp^ 


dxp^ 


r dxp^ 


(, M _L, ,0 „,iJ _ , 

Vj y^n,j ' ^ndkVk ^n,j) 


JJ 


(42) 


If one compares the definition of (see Equation [2^ to the expression in the bracket on the RHS of 
Equation |42l it yields 


ra ' o /n ra o /.n 




.IJ 


dp)l 




diPl 


dxj 


(43) 


As a result, one can replace the term /d'ip^)'ipl, in Equation 40, with the averaged expression 

derived in the above Equation Similarly to the derivation relevant to the translational equation 


of motion, one can replace the term {dwl ’’on the RHS of Equation 40, by the average value 
l/2{dwl'^/dxjn)ym for each facet. Equation 40 can be then rewritten as 




^ A(9wi\jj , dml'’ JJ 


I 


dxj 


Vi + 


dxj 


+ f E - MdmXA 


(44) 


18 


























Using = Xj + in the definition of wl’’ along with the identity equations d{ml'^nY)/dxj = 


dnif'^/dxjTiy and d{wf'^ny)/dxj = dwf"^/dxjny, Equation 


44 


can be written as 


PimiV 


V 


2Un 


EE ^{yj ^imk^mtk )j 

I Ti 


i^^ijk^jbk Pu^ijk^j^k) 


JL 

2Vo 


EE ^{Vj ^irnk^rn^k Vj ) 5 , 

I Ti 


(45) 


The last term on the RHS of Equation 45 is written considering the fact that 6° and i}° are equal for 


all particles inside the RVE. Furthermore, the first term on the RHS of Equation 45 can be expanded as 


{vYSimkXmtk’),j = SijkvYtY + eimkXm{yY^k’),ji wMch dyY/dxj = 0 is used. Therefore, Equation 
becomes 


45 


pL(r‘-:.)=^EE AsijkXj tj^ + EE A{xj Xj 771 

/ Xj I Xj 


^ + bl- Puvl 

^ Xi 


(46) 


The last term on the RHS of Equation is the moment of the translational equilibrium equation of 
the RVE (see Equation [3^ around the origin of the macroscopic global coordinate system; therefore, it is 
equal to zero. Comparing the first term on RHS of Equation [4^ with definition of macroscopic stress tensor 
of the RVE in Equation 38, one can replace it with The second term on the RHS of Equation 46 

is the divergence of the averaged moment stress tensor of the RVE. The macro-scale rotational equation 
of motion can be then written as follows 


f) ij^. 


(47) 


where 


= 

Pij 


^EE'4''K7s+‘"<?s) 


(48) 


I 


and the matrix Qf, is defined as Qf, = nYSjkix^eYi ■ X is the macroscopic moment stress tensor 


C4J ,,o 


19 











calculated using the results of RVE analysis, and Equation 1^ corresponds to the classical rotational 


equilibrium equation of Cosserat continua gsiiini. According to Equation ^ for macroscopic moment 
stress tensor and considering that + c^, one can conclude that consists of three terms: (1) the 

effect of the facet couple traction m; (2) the effect of the moment of the facet stress traction t around 
the particle node which the facet belongs to, and (3) the effect of the moment of the facet stress traction 
t, transferred to the particle node, around the centroid of the RVE. As result, the moment stress tensor 
is characterized by three length scales: (1) the facet size, associated to m; (2) the particle size or facet 
spacing; and (3) the size of the RVE. 


4 Numerical Results 

The homogenization theory formulated and discussed in the previous sections was implemented in the 
MARS computational software with the objective of upscaling the Lattice Discrete Particle Model 
(LDPM). LDPM, formulated, calibrated, and validated by Cusatis and coworkers [IZIIIH], is a meso-scale 
discrete model which simulates the mechanical interaction of concrete coarse aggregate pieces. LDPM has 
shown superior capabilities in modeling concrete behavior under dynamic loading [CT [52], Alkali Silica 
Reaction (ASR) deterioration [33] , as well as failure and fracture of fiber-reinforced concrete [331 [33] . 

The complete LDPM formulation is summarized in Appendix It is worth mentioning here that 
the LDPM computational units are polyhedral cells whose construction is anchored to the Delaunay 
triangulation of the simulated concrete aggregate pieces that are assumed to be spherical and size-graded 
according to the Fuller size distribution. In the LDPM formulation, each polyhedral cell represents one 
concrete spherical aggregate piece embedded in the surrounding mortar and the interfaces among the cells 
represent potential mortar cracks. Figure]^ shows a typical LDPM system of polyhedral cells and Figure 
Ht its periodic approximation. 

The generic RVE shown in Figure]^ is constructed as follows. Eight nodes are created at the vertexes 
of a cube (Figure]^). Then nodes are randomly placed on a RVE edge parallel to x axis, see node a in 
Figure]^. Then, these nodes are duplicated on the other three parallel edges along the x axis, see nodes 
b, c, and d in Figure]^. Similar procedure is carried out over the edges parallel to y and .2 axes. Next, the 
node generation on the RVE surfaces is performed by randomly placing nodes on a cube face with z axis 
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(a) (b) 


Figure 2: Polyhedral particle distribution in a LDPM prism: (a) generic LDPM system, (b) Periodic LDPM 
system. 



Figure 3: RVE generation procedure (a) Corner nodes (b) Edge nodes (c) Eace nodes 


as normal vector, see node e in Figure]^. The same nodes are then duplicated on the opposite RVE faces, 
see node / in Figure]^. Nodes on parallel cube faces with x and y axes as normal vectors are constructed 
with the same algorithm. Finally, nodes are placed inside the RVE based on the general LDPM procedure 
(see Appendix and relevant publications [IT] for details). 

As mentioned earlier in this paper the RVE analysis is conducted by imposing periodic boundary 
conditions. This is obtained by setting the displacements and rotations of the RVE vertexes to be zero 
and by imposing, through a master-slave constraint, that the periodic edge nodes and face nodes have the 
same rotations and displacements. 

The overall multiscale numerical procedure adopted in this paper can be summarized as follows. 

• The finite element method is employed to solve the macro-scale homogeneous problem in which 
external loads and essential BCs are applied incrementally. During each numerical step, strain 
increments = AuV — £jjfcA(p° and curvature increments A^^j = AcjT tensors are calculated at 
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each integration point based on the nodal displacement and rotation increments of the corresponding 
finite element. 


• The macroscopic strain and curvature increments are projected into the RVE facets through the 
proper projection operators: Ae^ = T*" {^lij + £jmn^>^imyn) = P^Anij. These projected 

strains and curvatures are imposed, upon sign change, as eigen-strains and eigen-curvatures, Ae° = 


Ae^ Ael = Ae^ - (-Ae^) and A-^^ = A'^jj, + Aipl = A'^jj, - (-A'^^) (See section [T^, to the RVE 
allowing the calculation of the fine-scale solution governed by the fine-scale constitutive equations. 


• Finally, the fine-scale facet tractions and moments are used to compute, through Equations and 
the macroscopic stresses, (j°-, and couple stresses, /rE for each Gauss point in the EE mesh. 
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4.1 Elastic RVE Analysis 

This section presents the analysis of the elastic macroscopic behavior of one LDPM RVE. The macro¬ 
scopic homogenized behavior is analyzed with reference to the classical constitutive equation for Cosserat 
elasticity, which, in non-dimensional variables, can be written as: 


Po^kk^ij T Pl^{ij) T P2'y[ij] j Pij Qo^kk^ij T T Q2^[ij] 

where foj = aij/(2fi + x) foj = L/ijj/[(2/i -|- are the normalized stress and couple tensors, L = 

characteristic size of the structure of interest, D = size of the RVE; Xij — liji >^ij — Lnij normalized strains 
and curvatures; po = A/(2/i-Fx), pi = l;p 2 = x/(2p-Fx); Qq = 7ri/[(2/i-Fx)-C>^], fo = (7r2-h7r3)/[(2/i-Fx)-C>^]; 

<12 = (7!'2 — 7r3)/[(2/i-|-x)T^^]; = kronecker delta; p, A, x, t<Ii t< 2 i and tts are the elastic constants; and the 

subscript parentheses and brackets represent extraction of the symmetric and antisymmetric, respectively, 
part of the tensors. 

In this section, eight different LDPM RVE sizes D= 15, 20, 25, 35, 50, 75, 100, and 150 mm are 
considered and 5 RVEs, characterized by different placement of the aggregate pieces, is studied for each 
case. It is worth mentioning that, in LDPM, different spherical aggregate placement inside the RVE 
yield to different RVE polyhedral particle configurations. The numerical calculations were performed by 
assuming the concrete mix design and model parameters reported in Appendix Figure shows the 
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homogenized values of po, P 2 , and of the normalized Young’s modulus defined as e = E/{2ii + x) = 
(3A + 2/i + x)/(2A + 2/i + x), as function of the RVE size normalized by the maximum spherical aggregate 
size, d = D/da- The error bars represent the scatter in the results obtained by simulating 5 different 
RVEs of the same size but with different realization of spherical aggregate positions inside the RVE. As 
one can see, the calculated values of the parameters tend to converge to a constant value as the size 
of the RVE increases and, at the same time, the results become independent of the spherical aggregate 
distribution inside the RVE. The value of p 2 is very close to zero for all RVE sizes and decreases rapidly 
with respect to the RVE size; this suggests that, for the analyzed fine-scale model, the homogenized stress 
tensor is symmetric. This result is due to the fact that in the LDPM formulation facet moments are 
zero, and this leads to facet traction distributions around each particle that have zero moment resultant 
around the particle node. In Figure the homogenized Poisson’s ratio is reported based on the equation 
u = A/(2A + 2ijl + x) and the calculated asymptotic value, 0.18, corresponds well with the value of 0.175 
calculated by exploiting the equivalence between particle models and microplane models nzi. Finally, 
Figure shows the homogenized parameters, go, and q 2 , as a function of the RVE size. These 
quantities also converge to an asymptotic value and become independent of the RVE spherical aggregate 
distribution for large enough value of D/da- By virtue of these results and by recalling the definitions of 
go, gi, and g 2 , it is interesting to note that the macroscopic Cosserat elastic parameters of the homogenized 
continuum depend quadratically on the RVE size. 

4.2 Nonlinear RVE Analysis 

In this section, the nonlinear response of the RVE is investigated under different strain and curvature 
loading conditions. Three different RVE sizes, D =25, 50, and 100 mm, and 7 different spherical aggregate 
placement inside the RVE are considered for each case. Typical polyhedral particle systems and geometry 
of each RVE size are shown in Figure]^ The nonlinear homogenized behavior of the RVE is studied under 
the effect of uniaxial strain tension and compression, hydrostatic compression, bending and torsional 
curvatures. In the following numerical examples, concrete mix design and model parameters are the same 
as the ones used in the elastic analysis. 
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Figure 4: Variation of elastic normalized effective material properties: (a) po, p 2 and normalized Yonng modnlus 
E. (b) u Poisson’s ratio, (c) qo, qi and q 2 , with respect to the ratio of RVE size to maximum spherical aggregate 
size. 

4.2.1 Nonlinear Analysis of RVE subject to components of the strain tensor 

Figure shows the homogenized stress-strain curves for different RVE sizes and polyhedral particle re¬ 
alizations relevant to RVEs subjected to uniaxial tensile strain. The results illustrate that the different 
polyhedral particle realizations do not affect the linear elastic and nonlinear pre-peak responses, but on the 
other hand, it clearly influences on the post-peak softening response. One can notice that the post-peak 
response of smaller RVE sizes is more scattered, while fine-scale randomness effect on the homogenized 
response diminishes for the larger RVEs 1231 EH. Therefore, one can conclude that the mesh realization is 
a more influential factor on the post-peak softening response of the RVEs of smaller sizes. Average of peak 
stress and strain values of different mesh realizations are calculated for each RVE size, and its variation 


with respect to the RVE size is plotted in Eigures [7a] and pb} As one can see these quantities as well as 
mesh realization effect decrease as the size of the RVE increases. 

Eurthermore, the average stress-strain curves of different polyhedral particle configurations for each 
RVE size are calculated and plotted in Figure As one can see clearly, increasing size of the RVE 
affects the post-peak behavior and increases the brittleness of the response. This is consistent with the 
well-known size effect associated to damage localization in quasi-brittle materials [53] • 

This phenomenon is depicted in Eigurej^ which shows damaged RVEs of different sizes at the end of the 
tensile loading process. The contour plots present meso-scale crack opening distributions corresponding 
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(a) 25 mm 


(b) 50 mm 


Figure 6: Macroscopic stress-strain curve for three different RVE size: 
tension 



(c) 100 mm 

25mm, 50mm, 100mm under uni-axial 


to macroscopic imposed uniaxial strain equal to 10 One can easily notice that the damaged area does 
not scale with the RVE size leading to the post peak size dependency on the RVE size. 


Evolution of damage for a 100 mm RVE is also shown in Figure [10] at five different macroscopic strain 
levels. Strain levels (1) and (2) are in pre-peak regime, in which damage is distributed throughout the RVE, 
which corresponds to the fact that homogenized response is not size dependent in the pre-peak regime. 
At strain level (3) which corresponds to the peak of the stress-strain curve, damage is still distributed 
over the RVE; However, as the material undergoes softening, damage localization initiates. Strain levels 
(4) and (5) are relevant to the softening branch of the response, in which damage localization is clearly 
visible. The size dependence of the homogenized softening RVE response leads to mesh-dependence of 


25 






















(a) (b) 

Figure 7: Variation of (a) average peak stress and (b) average peak strain, with respect to the RVE size. 



Figure 8: (a) Average tensile stress-strain curves for three different RVE sizes, (b) Coarse- and fine-scale strain 
energy density for different RVE sizes. 


the macroscopic response. This issue has been investigated by some authors [231 EH [59] with reference 
to continuum-based fine scale models. The complete analysis of this aspect with reference to the current 
LDPM-based homogenization scheme will be pursed in future work by the writers. 

Finally, in Figure [Sbl, the Hill-Mandel condition is verified by comparing the RVE strain energy density 
calculated through fine-scale and macroscopic quantities. 

Next, the nonlinear homogenized behavior of the RVE is studied under confined (uniaxial strain) and 
hydrostatic compression. For the confined compression test, a strain tensor with a longitudinal component 
up to -0.03 is considered, whereas for the hydrostatic compression case, all normal components of the 
strain tensor are set equal and with value up to -0.03. Figure [TT] shows the nonlinear response of RVEs of 
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Figure 9: Crack opening contour of damaged RVEs at tensile strain equal to 0.001 (left) 25 mm (middle) 50 mm 
(right) 100 mm 

different sizes and 7 different polyhedral particle configurations. 

In this case, due to the confinement, the stress-strain response is strain-hardening, and as one can 
see the different polyhedral particle realizations do not affect significantly the homogenized response in 
both the elastic and inelastic regime. In addition, the average of different mesh realization stress-strain 
responses is calculated and plotted for each RVE size in Figure |12a[ The nonlinear compressive response 
does not depend on the RVE size, which is consistent with the fact that plastic deformations are distributed 
through out the specimen, and strain localization does not take place. Finally, the Hill-Mandel condition 
is verified with reference to the confined compression test, and the fine- and coarse-scale strain energy 
density of different RVE sizes are plotted in Figure [T^ 

4.2.2 Nonlinear Analysis of RVE subject to components of the curvature tensor 

In this section, the nonlinear homogenized behavior of RVEs of 3 different sizes, 50, 75, and 100 mm and 
5 five different mesh configurations for each size, is studied under the effect of components of macroscopic 
curvature tensor. Bending and torsional behavior of the RVEs are investigated by applying macroscopic 
curvature tensors with the only non-zero components of ni 2 = 1 and kh = 1, respectively. Figure 
shows crack opening contour of damaged RVEs at the macroscopic curvature for Ki 2 = 0.5. The resulting 
crack pattern conforms with the fracture mode that one may expect from bending theories. Multiple crack 
lines are generated in the tensile strain domain, which is the top half of the RVEs, while half bottom 
part in under compression. As more strain is applied in compressive part, splitting cracks take place in 
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D = 100 mm 




Figure 10: Damage evolution of 100 mm RVE in tension test 


the latter region due to transverse tensile stress. Typical crack pattern of RVEs under torsion are plotted 
in Figure for Kh = 0.5. Crack opening contours show that the amount of damage close to the RVE 
center is negligible, while it increases as the facets are placed at further positions. This corresponds to the 
deformation mechanism and strain distribution in solids subject to torsion. 

Homogenized moment stress components /ii 2 and /in versus macroscopic curvature tensor components 
Hi 2 and Kii of RVEs of different sizes and polyhedral particle configurations are plotted in Figure [T^ One 
can see that effect of different polyhedral particle realizations on the homogenized response is negligible, 
which is due to the occurrence of distributed damage inside the RVE. Homogenized response of RVEs with 
different polyhedral particle realizations are averaged for each size and plotted in Figure p^Sb} It can be seen 
that the homogenized response consists of an initial elastic part and a hardening branch, which is related 
to the confinement due to the fact that all components of the strain tensor are zero, and the RVE cannot 
expand laterally. It is illustrated that, at any level of macroscopic curvature, magnitude of the moment 
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(a) D = 25 mm (b) D = 50 mm (c) D = 100 mm 

Figure 11: Volumetric stress-strain curve for three different RVE sizes under confined compression and hydro¬ 
static compression 


stress for larger RVE sizes is bigger compared to the smaller ones. Size dependency of moment stress was 


discussed in Section |4.1[ and it was shown that the elastic Cosserat coefficients are proportionai to the 
RVE size squared. In order to study size dependency in the nonlinear regime, the normalized quantities 
fiij = iiij/D and kij = kij x D are plotted in Eigure 


15c 


One can see that the normalized curves of three 


different RVE sizes are unique for both bending and torsion. This implies that the proportionality of the 
homogenized micropolar properties to the RVE size squared is still valid in the nonlinear regime. In Eigure 


15d, the Hill-Mandel condition is verified, and coarse- and fine-scale strain energy density are plotted for 


each RVE size for both the aforementioned cases. 

Einally, the existence of coupling effect characterized by the dependency of the homogenized stress 
tensor on the curvature tensor is investigated in both elastic and nonlinear regimes. The macroscopic 
curvature Ki 2 = 1 is applied on the RVEs, and the trace of the homogenized stress tensor is calculated and 


plotted in Eigure |15f| for different RVE sizes. One can observe that the trace of the stress tensor for the 
case of elastic RVE behavior is zero throughout the analysis. On the other hand, for the case of nonlinear 
behavior, it increases monotonically with the curvature. This implies that in elastic regime, stresses and 
strains are totally uncoupled from couple stresses and curvatures, whereas these quantities are strongly 
coupled in the nonlinear case. This aspect has been investigated very little in the literature where fully 
uncoupled behavior has been always postulated. 
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(a) (b) 


Figure 12: (a) Average compressive volumetric stress-strain curves for three different RVE sizes, (b) Coarse- 
and fine-scale strain energy density for different RVE sizes. 



Figure 13: Crack opening contour of damaged RVEs at bending curvature equal to 0.5 (a) 50 mm (b) 75 mm 
(c) 100 mm 


4.3 Tension Test on a Concrete Prism with Parallel Elastic Bars 


In this section, the behavior of a reinforced concrete prism under tension is studied in a full fine-scale 
simulation, and the obtained results are compared to the solution of the same problem through a two-scale 
homogenization algorithm, in which the concrete prism is modeled as a homogeneous continuum with a 
meso-scale material RVE assigned to every macroscopic integration point. Figure [T6a| shows the concrete 
prism and the two elastic bars attached to it, which are simulated by LDPM and solid finite elements, 
respectively. The same specimen in a two-scale homogenization problem is depicted in Figure |16b 


m 


which concrete prism is modeled by tetrahedral finite elements. Cross section of the concrete prism is 
100 mm X 100 mm, and its height is 500 mm. Two rigid loading plates are attached at the top and 
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Figure 14: Crack opening contour of damaged RVEs at torsional curvature equal to 0.5 (a) 50 mm (b) 75 mm 
(c) 100 mm 


the bottom of the whole specimen cross section to apply the boundary condition. Young modulus and 
Poisson’s ratio of the elastic bars are 28 GPa and 0.18, respectively. The same LDPM parameters used 
in the previous sections are adopted here. The specimen is pulled in the longitudinal direction up to a 
displacement equal to 0.7 mm. The RVE size is chosen to be 30 mm which approximately corresponds 
to the volume of each tetrahedral FE in the coarse mesh. This is done to mitigate the mesh-dependence 
due to the softening behavior of the RVE. The concrete prism and the elastic bars are connected through 
a master-slave algorithm. The numerical simulations of the coarse scale are performed by neglecting the 
couple stresses which are expected to be negligible for this particular application. 

The global force-displacement response of the full fine-scale and homogenization problems are plotted 
in Figure |16c[ Since concrete prism and elastic bars are tied and deform together during the loading 
process, distributed damage takes place through the whole specimen during the initial stages of the loading 
process, see Figure pT^. This damage state represents the linear elastic and the first hardening segment of 
the stress-strain response of the structure. The same damage state is captured through the homogenization 
procedure. Figure shows finite elements normal strain distribution along the loading direction through 
the specimen. One can see that the strain values are all in the same range, and no localization has 
occurred. The response of the full fine-scale and homogenization problems show excellent agreement in 
the elastic and the first hardening segment. As further deformation is applied on the structure, damage 
localizes in one section of the concrete bar and this causes a sudden drop in the global force-displacement 
curve. Subsequently, since the elastic bars and the concrete prism are forced to deform in parallel, the 
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(a) 


(b) 




(c) 



(d) (e) (f) 

Figure 15: (a) Homogenized couple stresses /in and /ii 2 versus curvature kh and K 12 for five different polyhedral 
particle configurations for each RVE size, (b) Average of the homogenized couple stress of different polyhedral 
particle configurations for each RVE size for kh and K 12 cases, (c) Scaled couple stress versus curvature curves, 
(d) Macro and Eine-Scale strain energy density evolution, (e) Scaled strain energy density evolution for the case 
ki 2 - (f) Trace of stress tensor due to elastic and nonlinear analysis of RVE under macroscopic K 12 . 


overall system can carry more load leading to a rehardening of the global response. Analysis of Figure 


16c shows that five damage localization events occur during the deformation process which corresponds 


to five sudden drops in the load-displacement curve. Crack pattern of the specimen is plotted after the 
formation of two, four, and five damage localization in Figure [T^, c, and d. It is interesting to show that 
the homogenization framework is able to generate the same damage distribution pattern. Figures [T7|F, g, 
and h show that two, four and five strain localization band appear in the specimen, which corresponds to 
the damage configuration obtained from full fine-scale problem. The global load-displacement curve of the 
homogenization problem also shows five sudden drops which conforms to the full fine-scale response, see 
Figure 16c[ The homogenized response captures well the displacement at which the first three localization 
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Figure 16: (a) Full LDPM concrete prism and attached elastic bars, (b) FE model of the concrete prism and 
attached elastic bars, (c) Force-Displacement curves obtained by homogenization and full fine-scale simulation. 


events occur, while it underestimates its value for later events. This is likely due to the relatively coarse 
mesh adopted at the macroscopic scale. 


5 Conclusions 

This paper presents the asymptotic expansion homogenization of fine-scale periodic discrete systems fea¬ 
turing independent translational and rotational degrees of freedom. Employing consistent asymptotic 
expansion of displacement and rotation fields, a rigorous analytical derivation was performed for elastic 
behavior, and it was extended to the nonlinear case upon making reasonable assumptions on the rigid 
body motions of a RVE. Based on this work, the following general conclusions can be drawn. 

• The equivalent homogenized continuum is of Cosserat-type characterized by nonsymmetric stress 
and couple tensors energetically conjugate to nonsymmetric strain and curvature tensors, respec¬ 
tively. The classical linear and rotational momentum balance equations can be derived from the 
homogenization of the fine-scale equilibrium equations. 

• The fine-scale kinematic quantities, namely facet strains and curvatures, are demonstrated to be 
related to the projection of the coarse-scale strains and curvatures into the local facet system of 
reference. This allows a straightforward implementation of the RVE problem into any computational 
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Figure 17: (Top row) Crack opening contour at different loading from full fine-scale simulation (Bottom row) 
Strain distribution contour at different loading steps from homogenization algorithm 


framework. 

• Similarly to previous research, the derived formula linking the fine-scale response to the coarse-scale 
stress tensor corresponds to the virial stress formulation commonly used for atomistic systems. 

• The derived formula linking the fine-scale response to the coarse-scale couple tensor is shown to 
consist of three terms with clear physical meaning. The first term is associated with the fine-scale 
couple tractions and it can be related to the facet size, which, in turn can be associated with the 
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size of weak spots in the material internal structure. The second term arises from the moment of the 
fine-scale stress tractions with respect to the particle node. As such, it depends on the size of the 
fine-scale particles and it can be related to the spacing or characteristic distance of the weak spots 
in the material internal structure. Finally, the third term is the effect of the moment of fine-scale 
stress tractions with respect to the center of the RVE and, consequently, it depends on the RVE size. 

The developed framework was then implemented in a computational software and applied to the upscal¬ 
ing of LDPM. Specific to this fine-scale model, the numerical results demonstrate the following interesting 
features of the equivalent homogenized continuum. 

• The macro-scale elastic parameters relating the stress tensor to the strain tensor become independent 
on RVE size and on the random position of the polyhedral particles inside the RVE for RVE sizes 
larger than about 5 times the maximum spherical aggregate size. On the contrary, the macro-scale 
parameters relevant to the relationship between curvature and couple tensors are shown to depend 
on the RVE size squared and they become independent on the random position of the polyhedral 
particles inside the RVE for RVE sizes larger than about 5 times the maximum spherical aggregate 
size. 

• The non-symmetric part of the macro-scale stress tensor is negligible since the relevant parameter is 
at least one order of magnitude smaller than the one governing the symmetric part. As a consequence, 
the linear and rotational momentum balance equations are decoupled. 

• In the elastic regime the stress-strain and couple-curvature constitutive equations are completely 
uncoupled. 

• In the non linear regime, for tensile loading and because the fine-scale behavior is strain-softening, the 
response is RVE-size-dependent. This is an expected result, although very often not acknowledged 
by most authors in the literature, associated with strain localization induced by softening. On the 
contrary, such dependence is not observed for compressive dominated loading conditions because the 
LDPM fine-scale behavior in compression is strain-hardening. 

• The coarse-scale couple-curvature constitutive equations scale with the square of the RVE size in the 
nonlinear range also but, contrarily to the elastic case, they show a strong coupling with the stress- 
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strain constitutive equations. Such coupling, never considered in the current literature of Cosserat 
media, will be studied in future work by the authors. 
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A Short Review of the Lattice Discrete Particle Model (LDPM) Geometrical Construction 
and Constitutive Equations 

LDPM model generation procedure and governing constitutive equations are explained in the following 
two sections. 

A.l LDPM model construction 

Concrete meso-scale structure is modeled by LDPM through the following steps: 

• Spherical aggregate generation is the first step which is carried out assuming that each aggregate 
piece can be approximated as a sphere. Under this assumption, the following spherical aggregate 
size distribution function proposed by Stroeven [56] is considered 


fid) = 


qdl 


[1 - (do/c^a)«]c?«+l 


(A.l) 


in which da and do are the maximum and minimum spherical aggregate size, respectively, and q 
is a material parameter. It can be shown jnS] that Equation A.l is associated with a sieve curve 


(percentage of spherical aggregate by weight retained by a sieve of characteristic size d) in the form 


= I f 


(A.2) 


where Uf = 3 — q. For Uf = 0.5 Equation A.2 corresponds to the classical Fuller curve which for its 
optimal packing properties, is extensively used in concrete technology. Considering concrete cement 
content c, water-to-cement ratio w/c, specimen volume, maximum da and minimum do spherical 
aggregate size along with the considered distribution function Equation |A.2[ the spherical aggregate 
system can be generated using a random number generator. 

• By using a try-and-error random procedure, spherical aggregate pieces are introduced into the con¬ 


crete volume from the largest to the smallest size. Figure [18a1 shows the spherical aggregate system 
generated for a typical dogbone specimen. 


42 









(a) (b) (c) 


Figure 18: (a) Spherical aggregate system for a typical dogbone specimen, (b) LDPM polyhedral particles for 
two adjacent spherical aggregate particle, (c) LDPM cell distretization for a typical dogbone specimen. 


• Delaunay tetrahedralization of the spherical aggregate piece centers is employed to define the inter¬ 
actions of the spherical aggregate system (Figure [l8b| ). 

• Finally, a three-dimensional domain tessellation anchored to the Delaunay tetrahedralization is car¬ 
ried out to create a system of polyhedral particles interacting through triangular facets, and a lattice 
system composed of the line segments connecting the spherical aggregate centers. Figure [T^ shows 
the final polyhedral particle discretization of a typical dogbone specimen. 


A.2 LDPM Kinematics 

The triangular facets forming the rigid polyhedral particles are assumed to be the potential material failure 
locations. Each facet is shared between two polyhedral particle and is characterized by a unit normal vector 
n and two tangential vectors m and 1. Accordingly, three strain components are defined on each triangular 
facet using Equations and which for LDPM gives 

_m^[ucl _F[ucl 

Cat —-; Cm —-; cl —- 

rp 

where [uc] is the displacement jump vector calculated at the facet centroid. One should consider that 
the LDPM constitutive equations explained in the next section are independent of facet curvatures. 
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A.3 LDPM constitutive equations 


This section reviews the specific constitutive equations governing the response of LDPM. First of all, it 
must be mentioned that LDPM assumes zero couple stresses at the meso-scale in both elastic and inelastic 
regime. This implies = 0 for a = A, M, L 

In the elastic regime, the normal and shear stresses are proportional to the corresponding strains: 
tiv = tM = Et^m] ti = Et^l, where Ej^ = Aq, Et = aEo, Eq = effective normal modulus, and 

a = shear-normal coupling parameter. Beyond the elastic regime, the vectorial constitutive relations are 
meant to reproduce three distinct sources of nonlinearity as described below. 

A.3.1 Fracture and cohesion due to tension and tension-shear 

For tensile loading (ejv > 0), fracturing and cohesive behavior due to tension and tension-shear are 
formulated through an effective strain, e = + Oi{e\j -\- e|), and stress, t = -|- {Im + thY/ot, 

which define the normal and shear stresses as = eN{t/e); Im = = aeL{t/e). The effec¬ 

tive stress t is incrementally elastic {i = E^e) and must satisfy the inequality 0 < t < crfct(e,a;) where 
= c’’o(<^) exp [—iLo(<^)(e — eo(<^))/c’'o(<^)], {x) = max{x,0}, and tan(a;) = j\foLeT = The 

post peak softening modulus is defined as Hq{uj) = Ht{2u>/7r)'^\ where rit is the softening exponent, Ht 
is the softening modulus in pure tension (a; = 7r/2) expressed as Ht = 2Eq/ [It/le — 1)', k — 2,EoGt/at] h 
is the length of the tetrahedron edge; and Gt is the mesoscale fracture energy. LDPM provides a smooth 
transition between pure tension and pure shear (a; = 0) with parabolic variation for strength given by 
ao{u}) = — sin(cc;) -|- y/sin^(a;) + 4:acos^{uj) /[2q; cos^(a;)], where Vgt — o's/o't is the ratio of shear 

strength to tensile strength. 

A.3.2 Compaction and pore collapse from compression 

For compressive loading (ejv < 0), the normal stress evolves incrementally elastically and is subjected to 
the inequality —<7bc{^DGv) < tv < 0 where ai,c is a strain-dependent boundary function of the volumetric 
strain, e^, and the deviatoric strain, en- The function expressing a^c models pore collapse for —ey < 
Gi = i^coGo = k^o-co/Eq, and it is formulated as a^c = (Tco + {-^v - ^co)Hc{rDv) where Hdrov) = 
-f^co/(l + ^c 2 {^Dv — ^ci)), o'cO IS the mesoscale compressive yield stress; and Kco, KcI, ^c 2 and 
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Hco are material parameters. Compaction and rehardening occur beyond pore collapse for —ey > e^i- In 
this case one has (Jbc = cfci{rDv) exp [(-ey - tci)Hc{rDv)/(^ciiTDv)] and (Jci{rDv) = (^co+{^ci-^ca)Hc{rDv)- 


A.3.3 Friction due to compression-shear 

The evolution of shear stresses simulate frictional behavior due to compression-shear. The incremen¬ 
tal shear stresses are computed as Im = ExieM — = Exi^L — e^), where = Xdi^/dtM, 

= Xdip/dtx, and A is the plastic multiplier with loading-unloading conditions < 0 and A > 0. The 
plastic potential is defined as Lp = \/t\j + — abs{tN), where the nonlinear frictional law for the shear 

strength is assumed to be abs = CTs + i/J^o — IJ'oo)o'm[^ — ^'^p{tN/<^No)] — fJ^ootN', is the transitional normal 
stress; /tq and Poo are the initial and final internal friction coefficients. 

Detailed description of model behavior in the nonlinear range can be found in Ref. m. 

A.4 Concrete Mix-Design and Model Parameters Used in the Numerical Simulations 

Minimum and maximum spherical aggregate size are do = 4 mm and da = 8 mm, respectively; cement 
content c = 612 kg/m^; water to cement ratio w/c = 0.4; aggregate to cement ratio a/c = 2.4; Fuller curve 
coefficient n/ = 0.42. 

The following LDPM parameters are used: Epf — 60 GPa, at — 3.45 MPa, Uco = 150 MPa, a — 0.25, 
fit = 0.4, It = 500 mm, Vgt = 2.6, Hco/Eq = 0.4, po = 0.4, = 0, kd = 1, kc 2 = 5, a no = 600 MPa, 

ol = Ex/E n — 0.25. 


B Asymptotic Expansion of Strains and Curvatures 


In order to obtain multiple scale definition of facet strain vector, one should first plug macroscopic Taylor 
series expansion of displacement and rotation of particle J around particle /, Eqs. and [T^ into facet 
strain definition, Equation In addition, equation x = r^y is used to change the length type variables to 
fine-scale quantities; = xy, rjcl = and r^c/ = c/. Equation writes 
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+ mjk ( + VO'imym + ^V^^imnym y'n ) 4 - mjkO]ci 
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2/3J 
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JJ 


(B.l) 


Spatial derivatives of displacement and rotation in equation above are partial derivative with respect 
to X. So, first and second order partial derivative of displacement and rotation asymptotic expansions, 
Eqs. [^and|^ with respect to x are as follows 


u. 




u, 




r]u, 


hj 


u. 


i^jk 


U, 


ijk 


mljt 


(B.2) 


0ij ~ rj + tp'lj + ojIj + Jjiplj Bijk ~ B + wljk 


(B.3) 


Using asymptotic expansion of displacement and rotation of a particle, Eqs. [^and[^ along with their 
macro-scale derivatives, Eqs. |B.2[ |B.3l and replacing them into Equation |B.1[ one obtains 
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Regrouping terms of the same order in above equation, one would get multiple scale definition of facet 
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In equation above, terms of order two and higher are neglected. Multiple scale definition of facet strain 
is derived, which consists of three classes of terms of 0{—l), O{0), and 0(1). Multiple scale definition of 
facet curvature vector will be obtained subsequently. Taylor series definition of rotation of particle J with 
respect to particle / in macro coordinate system. Equation should be inserted into definition of facet 
curvature. Equation]^ 


-1.^-1 


Xa = V r 


1 


oi + ri0YyY + Y^ Kikyryi - ^ 


JJ 


(B.6) 


Asymptotic expansion of rotation. Equation]^ along with its macroscopic first and second order deriva¬ 
tives, Equation |B.3[ are inserted into Equation |B.6| 


Xa = y J ^ 


v-^ujY + cpY + J + wY + + yJ^ylY + v^^mylY 


r7-a;°-^ /V''+ + Y^Y^ 


(B.7) 
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Collecting the terms of the same order and neglecting the ones of order more than zero, one can restate 
above equation as 
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Equation |B .81 is the multiple scale definition of facet curvature vector, which consists of terms of 0{—2), 
Oi-1), and 0(0). 


C Asymptotic Expansion of Facet Strain and Curvature using definition of rigid body 
motion of RVE 

Multiple scale definition of facet strain, Equation |B.51 and facet curvature, Equation |B.8[ can be rewritten 


regarding the definition of u°. Equation One can calculate first and second partial derivative of u° 
with respect to x as follows 




^ijk ^i,jk ' ^irnn^rn^jkyn 


(C.l) 


Using Eqs. C.l along with the fact that v°, A and A constant over the RVE: = v°, 

= A and A = A, one can revise Equation B.5 


€ry =r 




rj ^ 


+ EijkAYyi - vi) + EijuA^Y - Cfc)) 


+ I Uj Uj + SijkUJj SijkUJj Cf. 


+ vljVY + eijkAj{4 - 4) + ^imnl^^jy^Vn + Cfc^ 

+ v[4fjyY + Aliky^yk + Ain^nA^j^y^yiAi + AnkAY^nyY yy4 


+ £ijk I + Ajn^yY + AYr^yY) 4 - £ijkYY4 


^OLl 


(C.2) 


Using = y^ — y^ and = A — A in above equation along with y^^ + c'^ = y"^ , one would get 


48 

















erv = r 


^-1 


0/_lJ 




II-I I .,0 


,0.JJ 


,0 


V [Ui -Ui +eijkUJjC^-eijkUJjC^ + Vi jyj - EijkUJjVk + eijk^j^^y^Vk 
+ V ( uliv^ + £ijk^Y4 + ^ijk^^YmVY4 - ^ijkY/Y 


(C.3) 


+ T^YokY/yY + T^^ijk^lmnym yn yk + £ijk<^lmym 4 


0 


.0 ^.IJ^J 


JJ 


Multiple scale definition of facet curvature can also be revised by using = u;^ and = u;^ 

along with ~ can rewrite Equation Bi 


Xa = r ^ 


( ujY + 


+ y \ Vi +^i,jyj -^i +^i,jyj +7;^i,jkyj yk 


JJ 


(C.4) 


D Macroscopic Translational and Rotational Equilibruim Equations 

In order to derive macroscopic RVE translational equation of motion, one should consider the terms of 


0(1) in Equation 21 


Miu°' 


(D.l) 


Scaling back Equation D.l by multiplying both sides of the equation by y'^ and using the definition of 


t\ presented in Equation 29, one can get 


MLuf = n'£A^‘^\ + Vb« 


Ti 


dtl 


(D.2) 


where = YJ' Equation D.2 represents the 0(1) translational equilibrium equation for each particle 


inside the RVE. One can derive the RVE macroscopic translational equilibrium equation by summing up 
Equation |D.2| over all RVE particles and dividing by the RVE volume Vq 


-^MUv^^ + e- 

Vo 


imn^m Un^ 
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I Ti 
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e" + — W V^6° 


(D.3) 


49 
















In above equation is replaced by its definition, Equation 1^ Considering the fact that and uj'. 


, 0 / 




are equal for all RVE particles and the body force 6° is considered to be constant over the RVE, Equation 
ID .31 can be written as 




% E (y E ^4'!/^) = y E E E 


(D.4) 

wti V-o— / 

Second term on the left hand side of the Equation |D.4| is equal to zero considering the assumption 
that the local system of reference is the mass center of the particle system within the RVE; = 0. 

Einal form of the Equation |D .41 is presented in Equation |32| in Section [3!^ 

Macroscopic RVE rotational equation of motion can be derived by considering the terms of (9(1) in 
Equation!^ To have a consistent formulation for all particles and RVEs, one should consider the moment 
of all forces with respect to a fixed point in space, say the origin of a global coordinate system as shown 
in Eigure[^, which implies that the moment of Equation |D.l should be taken into account. Therefore, 
one can write 0{1) moment equilibrium equation of particle / as 


M^SijkY/ul^ + = X] ^ + QWai) + Y^^ijkYjbl 


(D,5) 


where is the position vector of particle / in the fine-scale global coordinate system Y = X/r7; 
Pa^cf = Y*" ^ moment of the facet traction with respect to the origin of the fine-scale global 

coordinate system, in which Y*^ = is the position vector of the contact point C between particles 

/ and J in the global coordinate system. Scaling back Equation |D.5| by multiplying both sides of the 


equation by 77^ and using the definition of p^, and presented in Equation 29, one can get 




+ ri-^Mj,Cjf = r,Y, A 

:fi 


( 1 

V 


dm^ 


d'ip\ 


^Pl]+V^e,,,Xhl 


(D.6) 


Equation D.6 represents the 0(1) rotational equilibrium equation for each particle inside the RVE. 
RVE macroscopic rotational equilibrium equation can be obtained by summing up Equation |D.6 over all 
RVE particles and dividing by the RVE volume Vq 
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T Y1 (vf+ekmnV ^ X] ^ = 

^ ^ Vo ^ 


IL 

K) 


(D.7) 


In above equation u^ is replaced by its definition, Equation 23. Considering equality of and uj\ 




for all RVE particles along with X| = Xj + xj, Equation 


D.7 


can be written as 


^ ^ T/ ^ ^ {^e^im + M^SijkSkmnXj^n) ^ 

yo j yo j 

X ^ ^ ^u^ijk^j^k ^ ^ ^ MuSijkSkmn^j^nV ~ 


T/ / ' Tr 

yo j yo j- 


(D.8) 
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Vo 


E E -4 (^4+i E ^ E 

Considering = 0 along with the equality of uf, Xj for all RVE particles, one can conclude 

that the third and the forth terms on the left hand side and the last term on the right hand side of the 
Equation |D.8| is equal to zero. Einal form of the Equation |D.81 is presented in Equation |40| in Section 
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